Reduced Operator Approximation for modelling open quantum systems 
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We present the Reduced Operator Approximation: a simple, physically transparent and compu- 
tationally efficient method of modelling open quantum systems. It employs the Heisenberg picture 
of the quantum dynamics, which allows us to focus on the system degrees of freedom in a nat- 
ural and easy way. We describe different variants of the method, low- and high-order (including 
the interaction operators), defining them for either general quantum harmonic oscillators baths or 
specialising them for independent baths with Lorentzian spectral densities. The wide range of the 
method applicability is presented on the example of systems coupled to different baths with differ- 
ent strengths, and compared with the exact pseudomode and the popular quantum state diffusion 
approach. Our results suggest that quantum coherence effects persist in open quantum systems for 
much longer times than previously thought. 
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I. INTRODUCTION 



The beginning of twentieth century launched a series 
of major paradigm shifts which heralded the era of mod- 
ern physics. It will perhaps be surprising to the mod- 
ern reader that in the advent of the revolutionary Ein- 
steinian theory of relativity, Maxwell and Boltzmann's 
kinetic theory and Planck's hypothesis of quanta, the sci- 
entific world was not convinced of the fact that matter is 
grainy and cannot be continuously divided infinitely 
The seed of doubt was planted by the renowned Scottish 
botanist, Robert Brown, who noticed in 1827 that pollen 
in water suspension which he examined under his micro- 
scope displayed a very rapid, irregular, zigzag motion. 
The mystery of the "vital force" driving the Brownian 
motions remained unsolved for nearly 80 years, evading 
the pincer of conventional physics. The answer came 
from Einstein and Smoluchowski, who showed how the 
behaviour of mechanical objects is driven by the statis- 
tical properties of thermal noise, postulating the exis- 
tence of molecules in the fluid and linking the diffusion 
strength of their motion to the friction acting on a body 
moving in the fluid @, Q • The explanation of Brown's 
experiments, being at the same time a major diversion 
from the "continuous" Newtonian dynamics forming the 
core of the contemporary physics, opened a whole new 
avenue of research into the behaviour of systems influ- 
enced with random noise, resulting in such fundamental 
discoveries as the fluctuation-dissipation theorem [3, [E[ . 
Since that time, dissipation has been shown to affect such 
key dynamical processes as electron transfer and trans- 
port, surface dynamics, quantum tunneling, control and 
nonadiabatic effects. More generally, scientists in many 
disciplines, from physics through biology to social sci- 
ences, have developed increasingly powerful methods of 
modelling open systems, which interact with their envi- 
ronment. 



In many nano-scale systems the noise influencing the 
dynamics arises from quantum fluctuations. Already in 
1928, when Nyquist proposed the fluctuation-dissipation 
theorem [J], the quantum fluctuations were treated dif- 
ferently than the classical ones: the energy fcgT from 
the classical equipartition law was replaced by the ther- 
mally averaged energy of a quantum harmonic oscilla- 
tor, a distinction becoming negligible at high temper- 
atures. This result has been followed by the devel- 
opment of the new branch of physics, the theory of 
open quantum systems ,6-8]. It has found applications 
in almost all areas of natural sciences [9j, from quan- 
tum optics [10| . through condensed matter physics [Til ], 
nanotechnology fl2| and spintronics through quan- 



15| . to biol- 
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turn information |14| . through chemistry 
ogy [l6[ or even stochastic gravity and inflationary cos- 
mology [l7| . Furthermore, it has implications for such 
fundamental problems as the quantum measurement the- 
ory UM and the emergence of classicality due to decoher- 
ence [19| . 

There exists a rich variety of methods of modelling 
open quantum systems, applicable to different physi- 
cal regimes and based on different approximation tech- 
niques (2014271 ] . In this paper we propose a new method, 
which describes finite-dimensional quantum systems in- 
teracting with non-Markovian quantum harmonic oscil- 
lator baths, from single modes to continuous spectra, 
as well as a wide range of interaction strengths, while 
having moderate computational requirements. The non- 
Markovianity is necessary to quantitatively analyse the 
properties of many physical systems encountered in the 
fields mentioned in the previous paragraph [TT|; 12814301 ] . 
The proposed method handles large or infinite baths and 
a wide range of interaction strengths, while having mod- 
erate computational requirements. It uses the Heisenberg 
picture, which makes it particularly easy to focus the at- 
tention on the system degrees of freedom while preserving 
the decoherence effects due to the coupling to the bath. 

In the following section we will remind shortly the the- 
oretical background of our research and the employed 
formalism fSecs. lITAl and ITXB|) . Next we will present the 
derivation of the Reduced Matrix Approximation ap- 
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proach ( Sec. Ill C|) and propose its two variants: low and 
high-order in the systems and bath operators. They 
will be optimised for typical cases of continuous and 
Lorentzian baths in Sec. Ill Dl In Sec. Mil we will present 
the results of our method and compare it to other known 
techniques of modelling open quantum systems, like the 
pseudomode method or the quantum state diffusion. Sec- 
tion llVI contains a short summary of our work. 



II. THEORETICAL APPROACH 

Most generally, an open quantum system is a subsys- 
tem of a larger, interacting quantum system, e.g. one 
of the photons in an EPR pair, an atom in a resonant 
cavity, a quantum dot interacting with phonons in the 
crystal or any real object "becoming classical" through 
scattering of a vast number of air molecules and photons 
on it. We consider the case of a finite-dimensional quan- 
tum system coupled to an infinite-dimensional quantum 
bath, composed of a possibly infinite number of modes. 
In such an asymmetrical setup it is natural to ignore the 
details of the bath dynamics and focus on the reduced 
density matrix of the system. In this chapter we derive 
this quantity using the proposed Reduced Operator Ap- 
proximation approach. 



A. Open quantum system 

We consider a quantum system represented in an N- 
dimensional Hilbert space H s spanned by basis states 
{|n)}. Its internal dynamics is described by the Hamil- 
tonian 

N 

H s ^ V-mntmn j 

m.n—1 

where t mn :— \m)(n\ are transition operators between the 
states \n) and |m) and V mn — V nm . In a more concise 
notation, H s is a trace of an N x N matrix product: 

H s = TrV T i, 

where t is a matrix of system operators, (t) rnn :— t mn . 
The Hermitian conjugate of a matrix of operators O is 
defined as (0*) mn := 0' nm . In the case of t, this leads to 
t tj since t mn — tmn- 

The system is coupled to a quantum bath composed of 
a collection of K independent harmonic oscillators living 
in an infinite-dimensional Hilbert space Jib, 

K 

fc=i 

where is the annihilation operator of the fc-th mode 
(h = 1). The coupling between the system and the bath 



is described by the operator 

K 

Hi = £(Tr&i)<4+h.c, (I) 

fc=i 

where (gk)mn = S mn gk n is an TV x N matrix describing 
the coupling of the fc-th bath mode with the system. The 
fact that each c/k is diagonal means that the bath does not 
induce transitions between system basis states. However, 
the matrix notation allows for easy generalisation of the 
model to include such bath-induced transitions. For any 
type of bath, we can define a spectral density function 

k 

The total Hamiltonian, generating the evolution of the 
system and bath in the Schrodinger picture, is given by 

H = H s + Hb + Hi = 

k k ^ 

TrV T t + y^a; fc a^afc + Tr^ {a\g k + a k glj t, 

k=l k=l 

where we have employed the fact that bath and system 
operators commute and v = t. 

The reduced matrix density of the system, p s (t), is 
defined as 

p s (t) := Tr b p(i) , 

where p(t) is the density matrix of the system and the 
bath as a whole, and the trace goes over the bath de- 
grees of freedom only. From the fact that operators t mn 
correspond to transitions between system basis states, it 
follows that p s (t) can also be computed from the formula 

p.(t)=Trp(t)i, (3) 

where the trace is applied to each operator in the p(t)t 
matrix of operators and taken over both system and bath 
degrees of freedom. The main task of the presented 
method is obtaining p s (t) without calculating p(t). 

B. Dynamics in the Heisenberg picture 

In the Heisenberg picture the wavefunction is time- 
independent, ^ = 'i'(O) (hence, the density matrix of a 
system in mixed state is time-independent as well), while 
an observable O (time-independent in the Schrodinger 
picture) satisfies 

j t O(t)=i[H,0(t)}, (4) 

where 0(t) := e lHt Oe~ lHt . From the last definition it 
follows that [Oi(t),0 2 (t)} = [0 1 ,0 2 ](t). 

We assume that at time t — the system and the 
bath — denoted by their initial reduced density matrices 
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p s and pb, respectively — are uncorrelated. Hence, p — 
p s ® pb and Eq. acquires the form 

p s (t) = Tr b i(t)p = Tr s [p s Tr b p fe £(t)] . (5) 

Let t(t) and a k (t) denote the Heisenberg-picture coun- 
terparts of t and dfc, respectively, with i(0) := t and 
ttfe(O) := Ofc. In order to derive their equations of mo- 
tions using Eq. (|4|) we calculate, for any A £ C NxN , 



[Tr(Ai{t)),t mn (t)} = 



N 

E 

m' ,?;/ — 1 

[A At)] 



Arn'n' \tn'm' (^) 5 ^mn (^)] 



(6) 



where we have used the identity (valid also in the 
Schrodinger picture) 



^mn(^)^m' n' (^) — trnn' (^)^rj 



(7) 



In the matrix-of-operators notation, [Tr(A£(t)), = 
[A, £] . The above identities can be applied to equation 
to obtain the evolution of t and a operators generated by 
the Hamiltonian ([2]). For the system we obtain 



K 

kt) = t[v T ,i(t)} + iY, (tikAtt)] + til hi) , (8) 



where s k (t) := t(t)ak(t) are system-bath interaction op- 
erators and, since v(i) — t(t) and bath operators com- 
mute with system operators, st(t) = i(t)al(i). For 
the bath, using the canonical commutation relations for 
bosonic creation/annihilation operators, we obtain 

flfc (*) = -iu k a k (t) - i Tr g k i (t) . (9) 



C. Reduced Operator Approximation 

1. General description 

The aim of the presented method is to model the evo- 
lution of the system, including its decoherence caused by 
the interaction with the bath. This information is con- 
tained in the reduced density matrix of the system p s (t) . 
As demonstrated in previous sections, Eqs ([3]) and ([5]), it 
can be obtained from the mean values of system operators 
t(t). Thus, to calculate p s (t) in the Heisenberg picture, 
one has to evolve t(t) in time. Additionally, since the 
evolution equation for t(t) ^ involves the bath operators 
dk(t), due to the system-bath interaction, it is necessary 
to evolve ak(t) as well. However, a numerical description 
of both groups of operators in the total system and bath 
basis is impossible, as %b is infinite-dimensional. 

According to Eq. (JSJ) , given the initial system state p s 
we only need to know the partial traces of system op- 
erators, TibPbt(t), to obtain p s (t). The corresponding 
partial traces of the bath operators, Tr b PbO-k(t), contain 



part of the information on how interaction correlated the 
bath and the system — if there was no such correlation, 
Tib Pb<ik (t) would be proportional to an identity opera- 
tor in % s . Thus, even after tracing out the bath de- 
grees of freedom, we can at least approximately capture 
the system-bath correlations arising from the interaction 
terms in Eqs. (|8|) and (HJ). This observation forms the 
basis of the proposed Reduced Operator Approximation 
(ROA). 

We represent both system and bath operators by N x N 
complex matrices in the system state basis (hence, t is 
represented by a matrix of matrices). Let M[0(£)] denote 
this reduced representation of an operator 0(t) with its 
elements defined as 



M[0{t)] mn :=Tr 6 p b (m\0(t)\n) 



(10) 



From the definition, M[0(t) t ] = M[0(i)]+. 

The evolution equations for the reduced representa- 
tions of system and bath operators are 

j t M[a k {t)] = -iio k M[a k (t)] - i Tr g k M[t\(t) (11) 



and 



-M[i(t)]=z[V T ,M[t(t)]] + 



K 

l^2([9 k ,M[sl(t)}) + [gl,M[s k (t)l 



fc=i 



(12) 



Since the system and the bath are correlated, 
M[s k (t)} ^ M[i(t)]M[a k (t)], which means that the above 
evolution equations are not complete. The simplest 
way to complete them is to approximate M[s k (t)] by 
the product of M[t(t)} and M[a k (t)]. However, again 
due to the system-bath coupling, M[i(t)]M[a k (t)] ^ 
M[a k (t)]M[i(t)\, even though [t(t),a k {t)] = 0. Thus, we 
need to specify a concrete ordering of the multiplied re- 
duced operators. We use the approximations of the form 

M[s k (t)} w 9M[t{t)]M[a k (t)] + (1 - 9)M[a k (t)}M[t(t)} , 
M[4(t)] ps eM[a k {t)]Hl[t(t)] + (1 - 6)M[i{t)]M[a k {t)]\ 

for 9 6 [0,1]. Numerical experiments have shown that 
simulations diverge for Q ^ \. Appendix [X] contains a 
simple theoretical explanation of this fact. Choosing 8 = 
\ , we arrive at the final form of evolution equation of the 
system operators in the reduced representation, 



-M[i(t)]=i[V T ,M[t(t)]} 



K 



J2\3k,{M\t(%M[a{(t)}}] 



k =i 

K 



(13) 



Y^[9l{M[t(t)},M[a k (t)]}}. 



k =i 



Equations (ITTj) and ([TBI employ reduced representations 
which are linear in the system or bath operators. Hence, 
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we will refer to them as the lower- order ROA. It is im- 
portant to note that this form retains the correlations 
between the system and the bath, due to the fact that 
bath operators are represented by their matrices in the 
system basis, M[a k (t)}. 

Additional information about the system-bath corre- 
lations is provided by the M[s k {t)\ matrix. Hence, it is 
beneficial to evolve it separately in addition to M[a k (t)] 
and M[t(t)}. For this purpose, we first derive the evolu- 
tion equation for s k (t), using Eqs. © and ©, 

j f h(t) = i(t)j t a k (t) + (jflt)) o*(*) 

K 



k'=i 

K 

i a l'(*)[fffc''*fc(*)] - iu k sk{t) - ii(t)g k 



fc'=i 



where we have used the fact that, due to the associativity 
of the operator product, §k>(t)ak{t) = i(t)a k >(t)a k (t) = 



a k '{t)s k (t) and 



sl(t)a k (t) 



aUt)S k (t). To de- 



rive the evolution equation for M[s k (t)], we have to 
solve a similar ordering problem as in the lower-order 
method. Differently than in the case of Eq. (H3J), in 
the evolution equation for M[s k (t)] we do not sym- 
metrise the reduced operator products, and obtain 
M[s k (t)a k ,(t)} « M[s k (t)]M[a k >(t)] and M[a£, (*)$*(*)] « 
M[a k ,(t)]M[s k (t)}. A simple argument justifying this 
approach can be found in Appendix [A] (we have also 
tested it for numerical stability). Furthermore, due 
to the aforementioned associativity of the operator 
product, M[s k (t)a k > (t)} can be approximated by either 
M[s k (t)]M[a k ,(t)], as above, or M[s k '(t)]M[a k (t)]. Sim- 
ilarly, we can choose between M[a k ,(t)]M[s k (t)} and 

M[§l,(t)]M[a k (t)] for M[a£, (*)«*(*)]• To exploit fully 
the information about the system-bath correlations con- 
tained in M[s k i(t)] matrices, we use an equally weighted 
average of the two approximations. In this way we obtain 
the evolution equation for the reduced representation of 
the interaction operator 



-M[s k {t)]=i[V T 1 M[s k {t)]] + 

K 



k' = l 



]T (\g k >,M[sUt)]} + \9l,M[s k .(t)]}) M[a k (t)\ 
J2M[a{ t (t)}l9k>,M[s k (t)}} ( 14 ) 



fc'=i 

K 



iJ2l9l,M[s k (t)}]M[a k ,(t)} 



k'=l 



- iui k s k (t) - ii(t)g k . 

Together with Eqs. (fTTj) and (fT2|) it defines the higher- 
order ROA. 



2. Reduced system density matrix 

From evolution equations (TT^I) or ([13)) we instantly see 
that, since the trace of every commutator is zero, 



N 



dt 



Tr M[t(t)}= ^2M[t mm (t)] =0. 



m— 1 



Inserting M[t(t)] instead of t(t) in Eq. ([5]), we obtain a 
trace-conserving expression for the reduced density ma- 
trix of the system, 



p.(t) =Tr p s M[t{t)} 



(15) 



However, if we use it to calculate p s (t), its positivity is 
not guaranteed. To fix this problem, we make use of the 
identity © to derive 

i(t)i(t) = Ni(t) 

and replace Eq. ([5]) with a different approximation 



Ps(t) 



N 



Tr p s M[t(t)] 



(16) 



Since t mm > (t)t m > n (t) = t mm >{t)(t nm >(t)y (and the same 
for the reduced representations), the above formula guar- 
antees that p s (t) is positive-semidefinite. On the other 
hand, the density matrix (TTB1) does not possess a con- 
served trace due to the fact that M [£(£)] 2 ^ M[i(t)i(t)}. 
Thus, we normalise the density matrix to obtain 



Psit) 



Tv Ps M[i{t)? 
Tr[Tv p s M[t(t)] 2 } 



D. Baths with continuous spectral densities 

In the limit of an infinite number of modes, spectral 
density J{uj) can be a continuous function. One way to 
handle this situation is to discretise J{uj) into a finite, 
but large number of modes. Assuming a constant mode 
frequency spacing Aw, we define a coupling constant for 
uj = kAuj to be 



9k 



AwJ(fcAw) 



(17) 



Taking the square root ensures proper normalisation as 
Auj — > 0. For a spectral density being a single Lorentzian 
peak, the method converges quite well already for K = 
100 modes per site and Aw ~ 7/IOO, where 7 is the half- 
width at half-maximum of the peak. 



1. Independent baths with continuous spectral densities 

When each system basis state is coupled to its own 
independent bath with a continuous spectral density, a 
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more sophisticated method is to describe these baths in 
terms of collective mode excitations caused by the cou- 
pling with the system. 

Let us consider bath operators multiplied by their 
phase factors, e Wkt a k (t), with dynamics given by 
j- t {e l ^ kt a k {t)) = -iTv g k i(t)e iuikt . Employing the fact 
that trace and integration commute, we obtain 



(t) = e~^ fet a fe (0) - % \ dse^ k{t ~ s ^ Tr g k t{s) 



which we insert into the definition of §k(t), 

h(t) = e- lLJkt i(t)a k (0) - ii(t) f dse'^-^ Trg k t(s) 

Jo 

= e- iUk H{t)a k {0) + u k {t) , 
and then the above formula into Eq. ([5]) , 

A" 

i(t) = i[v T ,i(t)} + iJ2 (&,«£(*)] + [glMt)} 



fc=l 



+ i 



1\ 

X; (e*"**4(0)[fffc^(*)] 



fc=l 



If each bath has an independent spectral density, each 
bath mode k is coupled to exactly one basis state n k , 

1-6- (gk)mn — 9knk^rn,nk^n,nk • Hence, 



N 



dt 



* r /■< 



E 



" l " fc(t ~ S) (|5fcm| 2 ^m(s) " |5fcn|^n«(s))^ - h.C. 



Sfcr, 



5^)e-^*a fe (0)+h.c. 



fe=i 



where we have used the fact that for independent densi- 
ties, g kn g~km 7^ only for n = m. In the limit of infinite 
number of modes, using Eq. (TIT)) we obtain 



A" 



lim Ve^^tl^^ff-i 

A — ^oo ^ — ' 



k=l 



whcr 



,(r) := / du)J m ((j)t 



is the bath correla- 



tion function. In this way we obtain a closed system 
of differential-integral equations for t mn (t), 



N 



^tmn(t) — ^ ^ ^ (Vm' nitm' 7i(t) Vnm'tmm'(t)) 



where 



a m (t) 



(18) 



H~ tmn(t) [\/ K m d m (t) K n d n {t) h.i 



A" 



'K m A->oo 
A 



lim y"g fcm a fc (t) 



— lim V^e^a^O) 

/ /^777 — tOO 



K poo 

and K m := lim V" |5fe m | 2 = / J m {u)<ko . 

fc=l 7-00 

Operators 5^(4) and a m (t) satisfy canonical commu- 
tation relations for bosons, 



[a m (t),a n {t)} = 0, 
[a m (t),a n (t)] = 



A 



lim y] g kn 9km = -5* 



v fc— 1 

They are pseudomode creation and annihilation opera- 
tors, creating or destroying collective excitations in a sin- 
gle bath (27j . Their dynamics is given by the equation 

d- M .. ^ u k e-^ kt a k {Q) q m (0) ^ 

a TO (t) = lim > Ofcm == 1 -^t mm (t) 



dt 



fe=i 



(19) 



Using the proposed method we have reduced signifi- 
cantly the number of bath operators, from K to N (for 
independent baths K > N, while in many cases K 3> N). 
However, numerical simulation of the differential-integral 
equation for the evolution of the reduced representation 
of d m (t) is difficult. In the next section we show that for 
a particular form of the spectral density function J m (w) 
one can get rid of the explicit time integration at the cost 
of a moderate increase of the number of simulated bath 
operators. 



2. Lorentzian spectral densities 

Continuous spectral densities composed of Lorentzian 
peaks, 



J n {uj) = X 



r, 



Inj 



IT (u - Ul nj ) 2 +7, 



(20) 



are especially popular due to their analytical tractability. 
In this section, we will optimise our method for this type 
of the system-bath coupling. The corresponding correla- 
tion function is a m {r) = r„ y e~ IWmj " r_7m; '' r '. Hence, 
for t - s > 0, 



dt 



iWmj(t— «)— 7mj (t— s) 



A continuous spectral density of the form (|20|) is con- 
structed from an infinite number of independent har- 
monic oscillator modes, with different modes contribut- 
ing to each Lorentzian peak. To derive the evolution 
equation for pseudomode bath operators we express them 
as sums of J2j Q-mj{t), where a m j(t) is constructed as 



mj 



lim 

A->oc 



gkmO>k\t) 



6 



and P? 1 C [1, . . . , K] is the set of indices of modes build- lution equation 
ing the j-th peak, UjPf = [1,...,K] and Pf n Pjf = 



Sjj'Pj 1 . This leads directly to [t mn (t), a m 'j< ft)] = 0, 
[a m j(t),a nj >(t)] = and [a\ nj {t),a n ,y(t)] = -8 mn 8jj>. 

Thus, fflt -(f) and a m j(t) are pseudomode creation 
and annihilation operators corresponding to individual 
Lorentzian peaks in bath spectral densities. 
Comparing a m j ft) with a m ft) leads to 



N 



^mj ft) 



lim V g^e-^akiO) 
r — ^nn * ■ 



Jo 



Differentiating over t gives 



ft) = lim V 



fcepf 



Ofc(0) 



(21) 



with the initial condition 



Omj(0) 



lim V" gkmak(0) ■ 



feepf 



By splitting a m ft) into a sum of a m:) - ft) , we have simpli- 
fied the differential-integral evolution equation ([19]) . In 
the reduced representation, 



dt 



M[a m j(t)] = ( -iu) m j-i/ m j)M[a m j{t)] + y/T mj M[t mm (t)] 

(22) 

with the initial condition M[a m j(0)] = 0. 

Since M[t mn {t)a m/j (t)} ^ M[t mn (t)]M[d m 'j ft)], 
we evolve separately the reduced representation of 
operator products s mnm 'j(t) := t mn {t)a m f j{t) = 

iT^j limK^ooI] fc 'ffto7s m „fe(t). Its adjoint equals 
s linm'jft) = *7im(*)2m'i(*)- Tne relevant commutators 
are [s mnm >j(t), a m »yft)] = and [s„ nm ^-(*), = 

Evolution equation of the system (|18p in the Lorentzian 
bath acquires the form 



d N 

^mn(^) — ^ ^ ^ (^m'm^m'nft) Vnm'tjnm' ft)) 
m' — 1 

+ V-Trnj S mnm j(t) — s nmm j(t) (23) 

+ E/ y/^nj s nmnj ft) — s mnnj ft) 
j 

Hence, the operators s rnnrn ij (t) themselves follow the evo- 



T.Smnm'j(t) — ^ ^ ^ (Vm"mSm"nm' j (t) Vnm" Smm"m' j ft)) 



ra" — 1 



+ E; [ s mnmj' ft) S nm-roj' ft) ^m'ift) 

V^™^ [ S Liij'W — Srnnnj'(t) d rn 'j(t) 

y 

+ t mn (t) lim > 



3 — iuifct 



fee p. 



K V 1 " 



-Ofe(0) 



(24) 



with the initial condition 

E/ 
3fcm' a fc(0) . 

Evolution equations for the reduced repretentations of 
the above system and interaction operators, respectively, 
are 

— M[t mn {t)] = i {Vm<mM[t m , n (t)] - V nm ,M[t mm ,{t)]) 
ra' — 1 

+ E VrVy [M[s mnmi ft)] - M[4 roroi ft)] 

i 

+ E v/r^J [^™ft)] - M[s mnnj (t) 

(25) 



and 



iV 



j.-^"[Smnm'j(^)] — 2 ^ ^ (Via' 1 rn^f\Sra !l nia' 



m" — l 



3 E V^^twWl - ^bi mm/ ft)])Af[a m ^(t)] + 
i' 

|M[8 mrem /j(t)] y^(\/r m j>M[a m j/ft)] - yT„j/M[a ni /ft)]) 
j' 

+ IE V^^^Lni'Wl - M[ S „-'ft)])M[a m ^ft)]- 

lE^Vrw^i^i'ft)] - \/iv M [«li'W])M[wi(i)] 



M[t mn (t)} 

(26) 

with initial condition M[s mnm ' 3 ■({))] — 0. We use the 
same operator order as in Subsec. HID II 

Higher-order Lorentzian ROA employs Eqs. (|22p . 
(|25|) and ([26]) . Lower-order Lorentzian ROA de- 
scribes the bath evolution with Eq. (|2"2"j) . To describe 
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the system evolution, we represent M[s mnm <j(t)] as 
^{M[t mn (t)],M[a m/j (t)}}, obtaining 



—M[t mn {t)] = i (V m > m M[t m > n (t)] - V nm ,M[t mm ,{t)]) 

m' — l 

+ ^{M[t mn (t)],J2[y/f~M[a mj {t)] 

3 

- ^jT~M[a nj {t)] -h.c.]}, 

(27) 

where the symmetrisation of the matrix product is justi- 
fied by the same arguments as in Sec. Ill C II 

The Lorentzian methods use much lower number of 
bath and interaction operators than the general ones, be- 
cause they model the bath excitations as collective pseu- 
domodes. Furthermore, thanks to the analytical integra- 
tion of the spectral density, they automatically include 
the tails of the spectral density, which are cut off by the 
discrete method. 



III. NUMERICAL EXAMPLES AND 
COMPARISON WITH OTHER METHODS 

In this section we present an example application of 
our method to the description of a molecular aggregate 
interacting with a non-Markovian quantum bath, and 
compare it with two other techniques: the pseudomode 
metho d |27j (PM) and non-Markovian quantum state dif- 
fusion (26| (QSD). 

The PM method replaces each Lorentzian peak in the 
spectral density by a pseudomode with a complex fre- 
quency, and models the dynamics of the original system 
and bath by simulating exactly, in the Schrodinger pic- 
ture, the system interacting with this pseudomode bath. 
As the reduced density matrix of the system p s (t) ob- 
tained in this way is exact, we use the PM method as 
our reference. The downside of the PM method is that, 
since it requires an exact simulation of a quantum many- 
body system, its computational requirements increase ex- 
ponentially with the number of system basis states and 
the number of bath spectral density peaks. The other 
method used for comparison, QSD, is an approximate 
one, which uses a Monte- Carlo simulation to calculate 
Pt(s). Its main advantage is the slow growth of the com- 
putational requirements with the system size. 

We model an exciton delocalised on a linear chain com- 
posed of N — 3 sites coupled by the nearest-neighbour 
potential, V mn = -(<5 m , n +i + £ TOj „_i). Each site in- 
teracts with a simple zero-temperature quantum bath 
with a unimodal Lorentzian spectral density J(w) = 
r77r~ 1 /((aj — w ) 2 + 7 2 ), where we set w = 1- We con- 
sider four cases: 




-1 -0.5 0.5 1 1.5 2 2.5 3 



bath mode frequency CO 

FIG. 1: Unimodal Lorentzian spectral densities used in sim- 
ulations. 
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The "wide" Lorentzian peaks correspond to fast decreas- 
ing bath correlation function, while the "narrow" ones 
indicate long correlation times. The coupling strength, 
"strong" or "weak" , determines the decoherence rate. 

We simulate the reduced density matrix of the system 
initially in the state ^ s — [1,0, 0] T , and compare the 
probabilities of finding the system in this state at later 
times, i.e. (p s (t))n. We use three variants of the ROA 
method: low-order ROA (Sec. Ill C II) . as well as low and 
high-order Lorentzian ROA (Sec. Ill D~2l) ; the high-order 
ROA has been skiped as the least efficient. The results 
are plotted in Figs. [2HH1 Comparison with the exact PM 
method shows that the best results are obtained for the 
Lorentzian variants of the method, which take into ac- 
count the whole range of the Lorentzian spectral density 
by analytical integration. At the same time QSD has the 
tendency to converge too rapidly to a steady state solu- 
tion, as compared to the exact PM method. The figure 
captions contain a detailed analysis of the results. 

For further comparison with the QSD method we cal- 
culate the transfer of the excitation on a ring aggregate 
in multimode Lorentzian bath (see Ref. [IH) in Fig. [5] 
The results obtained using the low-order Lorentzian ROA 
method are characterised by a much slower decoher- 
ence rate, as shown in the inset. This suggests that 
quantum coherence effects play a larger role in the dy- 
namics of excitons in open systems than predicted in 
Ref. [3l|. We attribute the faster decoherence observed 
in the QSD method to the zeroth order functional expan- 
sion (ZOFE) [32j required to make the method numer- 
ically feasible. It treats each path of the Monte Carlo 
simulation independently and thus increases artificially 
the amount of decoherence in the simulation causing the 
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FIG. 2: Comparison of the results of the ROA approach 
with the PM and QSD methods for bath A. The high-order 
Lorentzian ROA properly reproduced the amplitude and the 
phase of the probability density oscillations (the inset shows 
that it is superior to its low order variants) . The QSD method 
dampens the oscillations and does not reconstruct their phase. 



FIG. 4: Comparison of the results of the ROA approach 
with the PM and QSD methods for bath C. The high-order 
Lorentzian ROA precisely reconstructs the phase of the prob- 
ability density oscillations, while its low-order variants recon- 
struct the amplitude. The QSD method does not describe 
correctly the amplitude of the oscillations and loses the phase. 
The inset presents the comparison of the three ROA methods. 
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FIG. 3: Comparison of the results of the ROA approach with 
the PM and QSD methods for bath B. Although none of 
the methods satisfactorily reconstructs the PM results, the 
low-order Lorentzian ROA properly describes how the ampli- 
tude of fluctuations of the probability density changes in time. 
The inset presents a comparison of the three variants of ROA 
method; the high-order Lorentzian ROA diverges. 




FIG. 5: Comparison of the results of the ROA approach with 
the PM and QSD methods for bath D. The low-order ROA 
methods correctly describe the amplitude and phase of the 
probability density oscillations at shorter times and, together 
with the high-order Lorentzian ROA (which diminishes the os- 
cillation amplitudes), stabilise at the correct level. The QSD 
strongly diminishes the oscillations and fails to recover the 
correct steady state. 



system to converge too rapidly to a steady state solu- 
tion. Comparing the computational efficiency of ROA 
and QSD methods (where we have averaged the trans- 
fer over 1000 realizations of the stochastic noise [3l|), for 
the considered system our method is more than 10 times 
faster. 



IV. SUMMARY 

The presented Reduced Operator Approximation is a 
simple, physically transparent and computationally ef- 
ficient method of modelling open quantum systems. It 
employs the Heisenberg picture of the quantum dynam- 
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2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 
monomer monomer 



FIG. 6: Comparison of the results of the ROA Lorentzian 
method with the QSD method for a ring aggregate of 15 sites 
(initially only the site 8 is excited) in multimode Lorentzian 
bath (see Ref. [3ll ] for the description of the spectral density 
and Fig. 3 therein for the simulation parameters used in the 
plot). Inset: coherence defined as Tr p 2 calculated using both 
methods. 



ics, which allows us to focus on the system degrees of 
freedom in a natural and easy way. We have described 
different variants of the method: the low-order ROA, 
the high-order ROA (including the interaction operators) 
and their versions for Lorentzian baths. They have been 
applied to different systems (coupled to different baths 
with different strengths) . Comparison of our results with 
the exact pseudomode and the popular quantum state 
diffusion method favours the ROA approach, while the 
performance of ROA (especially the low-order case) is 
much higher than in the case of two other approachs. 
Furthermore, the ROA simulations of the exciton trans- 
fer on a ring aggregate suggest that quantum coherence 
effects play a larger role in the dynamics of open quantum 
systems than predicted by QSD [2a |. 

The method has been derived for the simplest case of 
linear coupling between the system and the bath. How- 
ever, its general version can be easily extended to higher- 
order couplings. Its another advantage over methods 
which do not use the Heisenberg picture (such as those 
used in our comparisons) is that a single simulation of 
the reduced system operators can be used to generate re- 



duced density matrices for an arbitrary choice of initial 
system state. 



Appendix A: Operator representation products 

Let us consider a Hamiltonian with N sites and N bath 
modes, V — and (gk)mn = gkS m ,kS n ,k- It is easy to sec 
that in this case the system operators are diagonal and 
t m m(t) = t mm (0). Hence, we can solve Eq. © explicitly, 
obtaining 



a m (t) 



e—™ 4 f a m (0) - ^(e**** - l)t mm (t) 



and, in the reduced representation, 

M[a m (t)} = - e-^MKnit)} 



(Al) 



(A2) 



This shows that even in this simple case, the M[t mn (t)] 
cannot in general commute with M[a k (t)]. 
Let us assume that at time t, the relation 

M[t mn {t)]M[t m , n ,(t)\ =8 nm ,M[t mn ,{t)] (A3) 

is preserved. Hence, 

M[t mn (t)]M[a k (t)} = -<U — (1 - e-^'^mnW] . 

'.In 



M[a k {t)]M\t mn {t)] = -6 mk ^(l - e 



" lulmt )M[t mn {t)] . 
Applying the above results to Eq. (TP3")) . we obtain 

(A4) 



4fM[t mn (t)] = i(l - S mn )[d(a m (t) - a n (t)) 



+ (1 - 9)(a rn (t) ~ a n (t))]M[t mn (t)} 



-(1 — e *">"*). It easy to see that for 



where a m (t) 

n 7^ m', 4f.M[t mn (t)}M[t m i n i(t)} — 0, consistently with 
Eq. (|A3I) . For n = m', we have 



iM[t mn (t)]M[t nn >(t)] = i(l - 5 mn )[6(a m (t) - a n {t)) 



dl 



+ (1 - 6){a m (t) - a n (t))]M[W(*)] 



+ i{l-5 nn ,)[6{a n {t)-a n ,{t)) 



+ (1 - 0)(a n (t) - a n ,{t))]M[t mn >{t)] , 
while in the case of m = n, 



f t M[t mn {t)]M[t nn ,{t)] = i(l - S mn ,)[0(a m (t) - a„»(t)) 



+ (1 - 6)(a m (t) - a n ,{t))]M\t mnl (t)] = -M[t mn> (t)] , 

again consistently with Eq. (|A3I) . Analogous result is 
obtained for n — n' . However, in the case of m ^ n and 
ri 7^ n' , we have 



£ t M[t mn (t)}M[t nn ,{t)} = i[6(a rn (t) - a n (t))- 



(1 - 9){a m {t) - a n {t))]M[t mn ,{t)]+i[e{a n {t) - a n ,{t)) 



(1 - 9)(a n (t) - a n ,(t))]M[t mn , (t)] 
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It is only for 8 — 1 /2 that the terms with a n (t) cancel out 
and the right side is equal to ■^M[t mn <(t)]. Therefore, 
for 9 = 1/2 evolution equations for M[t mn (t)] preserve 
product identities for t mn (t). Consequently, we obtain 

f t M[t mn {t)] = i [Re (a m (t) - a„(t))] M[t mn (t)} , 

leading to the solution 

M[t mn (t)} = M[t mn (0)] exp(i / Re (a m (s) - a n (s)))ds. 

Jo 

For operators s mnk {t), we obtain from Eq. (IA1 [) that 
M[s mnk (t)} = -^6 nk (l - e- iu ^)M[t mn (t)} , (A5) 

assuming the bath was in the ground state at t = 0. 
For the simple model considered in this Appendix, their 
Heisenberg equations are 

+ iigmalnit) - 9nai,{t))s mn k{t) (A6) 

It is easy to notice that s mnk (t)a k >(t) = s mnk , (t)a k (t) 

and a\,{t)s mnk {t) = 4mfc'(*H(*)- From 
Eqs. (| A2|) and (|A5I) we see that this property 
is preserved in the reduced representation and 

M\s mnk {t)}M[a kl {t)] ~ 5 

nk^nk' 

M[t mn (t)} while 

M[al,{t)}M[s mnk (t)} ~ S 

nk^mk' 

M[t mn (t)}. Thus, 
when choosing the reduced representation of Eq. (|A6|) it 



is irrelevant whether the index k or k' is attached to the 
s operator on the right side. For simplicity, we keep the 
k index with the s operators, and use the representations 

M[s mnk (t)a k ,(t)} = 6M[s mnk (t)]M[a k >(t)] 

+ {l-d)M[a k ,]M[s mnk {t)], 

M[al(t)s mnk (t)} = 0M[a k , M[s mnk (t)] 

+ (l-6) (M[s mnk (t)]M[a k ^ - M[t mn (t)]6 k ,k>) , 

where 9 € [0, 1]. Applying this to Eq. (|A6I) . we obtain 

f t M[s mnk (t)} = i6M[s mnk {t)]{j^M[a m {t)]-g^M[a n {t)\) 
+ i9(g m M[al(t)} - g n M[<f n (t)])M[s mnk (t)] 
+ i(l - 9)(g^M[a m (t)\ - g^M[a n {t)])M{s mnk {t)\ 
+ i(l - 8)M[s mnk (t)}(g m M[al(t)} - g n M [oj, (*)]) 

- i(l - 9){g rn 5 rntk - g n 5n, k )M[t mn (t)] 

- iuj k M[s mnk (t)} - iS nk g n M[t mn (t)} . 

Since for m = n the ^-dependent parts are zero, we 
have to consider the case m ^ n to determine the best 
choice of 9. Additionally, we assume that k = m. 
From Eq. (|A3|) we know that M[s mnm (t)] = 0. Thus, 
consistency requires that -^M{s mnm (t)] calculated us- 
ing the above expression should evaluate to zero. Us- 
ing Eq. (IA2p and the above assumptions, we transform 
it to f t M[s mnm {t)} = - 9)g m M[t mn (t)]. Thus, con- 
sistency requires that we choose 9 = 1 to calculate the 
derivatives of M[s mnk (t)]. 
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